/**
  SUNNY project, Anyang Normal University, IMP-CAS
  \class TAZBM
  \brief First try of toy shell model using ZBM potential. Ref.: PRL.21.39
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2022/11/28
  \date Last modified: 2022/11/28 by SUN Yazhou
  \copyright 2020-2021 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <cstring>
#include <algorithm>
#include <iostream>
#include "TAMBSDManager.h"
#include "TASPStateManager.h"
#include "TAMBSDList.h"
#include "TAZBM.h"
#include "TAMBSD.h"
#include "TAException.h"
#include "TASPState.h"
#include "TACGCoe.h"
#include "TAMath.h"
#include "TABit.h"
#include "TAHF.h"

using std::min;
using std::max;
using std::swap;

static auto mi = TAMath::minus; // mi(n) = (-)^n
static auto cg = TACGCoe::CG;

TAZBM::TAZBM(const string &cfgFile) : TAHamiltonian(cfgFile){
  memset(fTBME, 0, sizeof(fTBME));
  // fTBME[6][6][6][2]: [ll][ll][J][T]
  // [6][6]: dd,ss,pp,ds,dp,sp; [6]: 0,1,2,3,4,5; [2]: 0,1 //
  // [dd][dd]
  fTBME[0][0][0][1] = -2.81;
  fTBME[0][0][1][0] = -1.30;
  fTBME[0][0][2][1] = -0.98;
  fTBME[0][0][3][0] = -1.02;
  fTBME[0][0][4][1] = +0.12;
  fTBME[0][0][5][0] = -3.86;
  // [ss][ss]
  fTBME[1][1][0][1] = -2.28;
  fTBME[1][1][1][0] = -4.04;
  // [pp][pp]
  fTBME[2][2][0][1] = -2.37;
  fTBME[2][2][1][0] = -4.55;
  // [dd][ss]
  fTBME[0][1][0][1] = -1.20;
  fTBME[0][1][1][0] = -0.93;
  // [dd][pp]
  fTBME[0][2][0][1] = +3.37;
  fTBME[0][2][1][0] = -1.50;
  // [ss][pp]
  fTBME[1][2][0][1] = +0.73;
  fTBME[1][2][1][0] = -0.50;
  // [dd][ds]
  fTBME[0][3][2][1] = -0.84;
  fTBME[0][3][3][0] = -1.69;
  // [ds][ds]
  fTBME[3][3][2][0] = -0.80;
  fTBME[3][3][2][1] = -1.15;
  fTBME[3][3][3][0] = -3.90;
  fTBME[3][3][3][1] = +0.24;
  // [dp][dp]
  fTBME[4][4][2][0] = -4.65;
  fTBME[4][4][2][1] = +0.67;
  fTBME[4][4][3][0] = -2.71;
  fTBME[4][4][3][1] = -0.95;
  // [sp][sp]
  fTBME[5][5][0][0] = -3.17;
  fTBME[5][5][0][1] = +0.35;
  fTBME[5][5][1][0] = -3.01;
  fTBME[5][5][1][1] = +0.47;
} // end constructor

TAZBM::~TAZBM(){}

// retrieve two-body matrix elements //
double TAZBM::TBME(int p, int q, int r, int s) const{
  static TASPStateManager &sm = *TASPStateManager::Instance();
  TASPState *sp = sm[p], *sq = sm[q], *sr = sm[r], *ss = sm[s];
  // 0. collect single-particle states information //
  const int lp = sp->l(), lq = sq->l(), lr = sr->l(), ls = ss->l();
  if(mi(lp+lq) != mi(lr+ls)) return 0.; // inconsistent parity
  const double mp = sp->mj(), mq = sq->mj(), mr = sr->mj(), ms = ss->mj();
  if(mp+mq != mr+ms) TAException::Error("TAZBM",
    "TMBE: Inconsisten M occurred.");
  const double tp = sp->tz(), tq = sq->tz(), tr = sr->tz(), ts = ss->tz();
  if(tp+tq != tr+ts) TAException::Error("TAZBM",
    "TMBE: Inconsistent Tz occurred.");
  const double np = sp->n(), nq = sq->n(), nr = sr->n(), ns = ss->n();
  const double jp = sp->j(), jq = sq->j(), jr = sr->j(), js = ss->j();
  const double M = mp+mq, t = tp+tq; double me = 0.;
  // alike: whether the two are at the same sp orbit
  const bool alikepq = (np == nq && lp == lq && jp == jq); // p(nlj) == q(nlj)
  const bool alikers = (nr == ns && lr == ls && jr == js); // r(nlj) == s(nlj)

  // 1. Get the allowed range of J //
  // [Jmin,Jmax] is the smaller one of [|jp-jq|,jp+jq] and [|jr-js|,jr+js]
  const double Jmax = min(jp+jq, jr+js), Jmin = max(fabs(jp-jq), fabs(jr-js));

  // 2. Sum of coupled TBMEs weighted by C-G coupling coefficients //
  for(double J = Jmin; J <= Jmax; J++) for(double T = 0.; T <= 1.; T++){
    // J+T should be odd to enforce antisymmetry if (p,q) or (r,s) are in the
    // same sp orbit (have the same nlj). NOTE that this requirement does not
    // apply to the situation where p and q are different because antisymmetry
    // could always be taken care of by rearranging the radial parts of the
    // wavefunctions, namely, R_{nlj}(p) and R_{nlj}(q)
    if((alikepq || alikers) && !(int(J+T) & 1)) continue;
    if(J < fabs(M) || T < fabs(t)) continue; // cgs are zero for these (J, T)
    const double cgJpq = cg(jp,mp,jq,mq,J,M);
    const double cgJrs = cg(jr,mr,js,ms,J,M);
    const double cgTpq = cg(0.5,tp,0.5,tq,T,t);
    const double cgTrs = cg(0.5,tr,0.5,ts,T,t);
    const double tbme = TBME(lp,lq,lr,ls,J,T);
    const double tmp = cgJpq*cgJrs*cgTpq*cgTrs*tbme;
    me += tmp;
    // me += cg(jp,mp,jq,mq,J,M)*cg(jr,mr,js,ms,J,M)*
    //   cg(0.5,tp,0.5,tq,T,t)*cg(0.5,tr,0.5,ts,T,t)*
    //   TBME(lp,lq,lr,ls,J,T);
  } // end for over J and T
  /// FIXME: May not be necessary ///
  // if(alikepq && alikers) me *= 2.; // DEBUG
  // else if(alikepq || alikers) me *= TAMath::Sqrt2(); // DEBUG
  /// oooOOOoooOOOooo ///
  return me;
} // end member function TBME

// Two-body matrix element
// fTBME[6][6][6][2]: [ll][ll][J][T]
// [6][6]: dd,ss,pp,ds,dp,sp; [6]: 0,1,2,3,4,5; [2]: 0,1
inline int hash(int l1, int l2){ // hash table: l1,l2,l3,l4->TBME
  if(2 == l1 && 2 == l2) return 0; // dd
  if(0 == l1 && 0 == l2) return 1; // ss
  if(1 == l1 && 1 == l2) return 2; // pp
  if(2 == l1 && 0 == l2) return 3; // ds
  if(2 == l1 && 1 == l2) return 4; // dp
  if(0 == l1 && 1 == l2) return 5; // sp

  TAException::Warn("TAZBM", "hash: illegal l1 and l2");
  return -1;
} // end inlinle function hash
double TAZBM::TBME(int l1, int l2, int l3, int l4, int J, int T) const{
  const int p = hash(l1, l2), q = hash(l3, l4);
  // if(l1 < l2) TAException::Warn("TAZBM",
    // "TBME: Odd, l1 < l2, l1: %d, l2: %d", l1, l2);
  if(!fTBME[p][q][J][T]) TAException::Warn("TAZBM", "TBME: Encountered 0.");
  return fTBME[p][q][J][T];
} // end member function TBME

double TAZBM::MatrixElement1N(int rr, int cc) const{
  // if rr and cc have more than 1 sps different, return 0. //
  TABit br((*fMBSDListM)[rr] ^ (*fMBSDListM)[cc]);
  if(br.count() > 2) return 0.;
  // the \sum_n{delta*n*N_n} term //
  if(rr != cc) return 0.;

  TAMBSD m((*fMBSDListM)[rr]); // DEBUG
  return m.e(); // sum of energies of all the sps's
} // end member function MatrixElement1N
double TAZBM::MatrixElement2N(int rr, int cc) const{
  // if rr and cc have more than 2 sps's different, return 0. //
  TABit br((*fMBSDListM)[rr] ^ (*fMBSDListM)[cc]);
  if(br.count() > 4) return 0.;

  double phase = 0.;
  if(br.count() == 4){ // two sps's different
    TABit brr((*fMBSDListM)[rr]), bcc((*fMBSDListM)[cc]); // the two sps's
    // find the two sps's by their subscripts //
    int a[2]{}, b[2]{}, cnta = 0, cntb = 0;
    for(int i = 0; i < fNSPState; i++)
      if(brr[i] ^ bcc[i]){ if(brr[i]) a[cnta++] = i; else b[cntb++] = i; }
    int p = a[0], q = a[1], r = b[0], s = b[1];
    // XXX ds,dp,sp (descending) are positive order, but sd,pd,ps (ascending)
    // is positive in default in this code. So additional phase is needed XXX
    if(p < q) swap(p, q);
    if(r < s) swap(r, s);
    // Integral(rr, p, q, s, r, cc): <rr|a+_p*a+_q * a_s*a_r|cc>
    // i.e. the order does not change in Integral(..), as in the operators
    phase = fMBSDListM->Integral(rr, p, q, s, r, cc);
    if(!phase) return 0.;
    return phase*TBME(p, q, r, s);
  } // end if
  if(br.count() == 2){
    // cout << "--------------  BINGO  ------------------" << endl; // DEBUG
    double me = 0.;
    TABit brr((*fMBSDListM)[rr]), bcc((*fMBSDListM)[cc]);
    // cout << "rr: " << rr << ": "; brr.Print(); // DEBUG
    // cout << "cc: " << cc << ": "; bcc.Print(); // DEBUG
    int p0 = -1, r0 = -1; // the two sp states that are different in rr and cc
    for(int i = 0; i < fNSPState; i++)
      if(brr[i] ^ bcc[i]){ if(brr[i]) p0 = i; else r0 = i; }
    if(-1 == p0 && -1 == r0)
      TAException::Error("TAZBM", "MatrixElement2N: shit happens.");
    for(int q0 = 0; q0 < fNSPState; q0++){
      if(q0 == p0 || q0 == r0) continue; // Pauli exclusion
      // use delegates of p,q,r,s in case of variable swaps
      int p = p0, q = q0, r = r0, s = q0;
      // swap (p,q) and/or (r,s) if they're not of the 'normal' order
      // 'normal' is what is given in ZBM TBMEs
      if(p < q) swap(p, q);
      if(r < s) swap(r, s);
      phase = fMBSDListM->Integral(rr, p, q, s, r, cc);
      if(!phase) continue;
      // cout << "DING~" << endl; // DEBUG
      const double mm = TBME(p, q, r, s);
      if(!mm) continue;
      me += phase*mm;
      // cout << "p: " << p << " q: " << q; // DEBUG
      // cout << " r: " << r << " s: " << s << endl; // DEBUG
      // cout << "phase: " << phase << " mm: " << mm << endl; // DEBUG
      // cout << "me = " << me << endl; // DEBUG
      // getchar(); // DEBUG
    } // end for
    return me;
  } // end if

  // br.count() == 0 //
  double me = 0.;
  for(int p = 1; p < fNSPState; p++) for(int q = 0; q < p; q++){
    // Integral(rr, p, q, s, r, cc): <rr|a+_p*a+_q * a_s*a_r|cc>
    // i.e. the order does not change in Integral(..), as in the operators //
    phase = fMBSDListM->Integral(rr, p, q, q, p, cc);
    if(!phase) continue;
    me += phase*TBME(p, q, p, q);
  } // end for over creation operators a+_p
  return me;
} // end member function MatrixElement2N

double TAZBM::MatrixElement3N(int rr, int cc) const{
  return 0.;

  // if rr and cc have more than 3 sps's different, return 0. //
  TABit br((*fMBSDListM)[rr] ^ (*fMBSDListM)[cc]);
  if(br.count() > 6) return 0.;

  return 0.;
} // end member function MatrixElement3N
